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Abstract 

A new computer program has been developed called ASP3D (Advanced Small 
Perturbation - 3D), which solves the small perturbation potential flow equation in an 
advanced form including mass-consistent surface and trailing wake boundary conditions, 
and entropy, vorticity, and viscous effects. The purpose of the program is for unsteady 
aerodynamic and aeroelastic analyses, especially in the nonlinear transonic flight regime. 
The program exploits the simplicity of stationary Cartesian meshes with the movement or 
deformation of the configuration under consideration incorporated into the solution 
algorithm through a planar surface boundary condition. The ASP3D code is the result of 
a decade of developmental work on improvements to the small perturbation formulation, 
performed while the author was employed as a Senior Research Scientist in the 
Configuration Aerodynamics Branch at the NASA Langley Research Center. The 
ASP3D code is a significant improvement to the state-of-the-art for transonic aeroelastic 
analyses over the CAP-TSD code (Computational Aeroelasticity Program - Transonic 
Small Disturbance), which was developed principally by the author in the mid-1980s. 
The paper presents unsteady aerodynamic and aeroelastic applications of ASP3D to 
assess the time dependent capability and demonstrate various features of the code. The 
cases considered include: (1) the NACA 0012 airfoil undergoing a forced harmonic 
pitching motion about the quarter chord at transonic conditions, (2) a thickening-thinning 
parabolic arc airfoil at a transonic Mach number, (3) an investigation of wave 
propagation characteristics with emphasis on high wave number applications, (4) 
aeroelastic transient calculations for the NACA 0012 airfoil at values of the dynamic 
pressure below, near, and above the flutter value, and (5) the F-5 fighter wing undergoing 
a rigid pitching motion about the wing root midchord axis. The results compare well 
with alternative methods and experimental data, and thus demonstrate the efficiency and 
utility of the ASP3D code for various unsteady aerodynamic and aeroelastic applications. 

Introduction 

Over the years research has been conducted by several organizations including the 
NASA Langley Research Center to assess the general applicability and accuracy of 
transonic small perturbation (TSP) computer programs in numerous unsteady 
aerodynamic and aeroelastic applications. 1 " 4 A notable example of such a computer 
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program is the CAP-TSD code 5 developed by a team of researchers at NASA Langley 
Research Center in the mid-1980s. The CAP-TSD code solves the transonic small 
disturbance (TSD) equation with a time-accurate approximate factorization (AF1) 
algorithm ' and couples the flow equations with the structural equations of motion for 
simultaneous time integration for aeroelastic analyses. The code exploits the simplicity 
of stationary Cartesian meshes with the movement or deformation of the configuration 
under consideration incorporated into the solution algorithm through a planar surface 
boundary condition. However, applications with the code have revealed various 

o 

inaccuracies inherent in the underlying small perturbation theory including accuracy 
limitations near the leading edge and wing tip regions, 1 shock waves that are inaccurately 
captured in terms of strength or location, 4 and convergence difficulties attributable to 
certain computational and mathematical complications such as non-uniqueness of the 
potential flow equation. 9,10 Consequently, the author conducted research over the last 
decade to examine the applicability of the small perturbation concept in general and the 
accuracy of the CAP-TSD code in specific. The author is the principal developer 6,7 ’ 11 ’ 12 
of the CAP-TSD code and is therefore in a unique position to make a critical assessment 
of the methodology contained therein. The objective of the effort was to identify the 
sources of the inaccuracies and determine if improvements could be made to alleviate or 
eliminate them. 

Subsequently, a new advanced small perturbation (ASP) potential flow theory 13 was 
developed by methodically determining the essential elements required to produce 
accurate solutions with the small perturbation approach on a Cartesian grid. The ASP 
theory involves a higher-order streamwise flux in the governing equation, mass- 
consistent surface boundary conditions for flow tangency, mass-conserving entropy and 
vorticity effects including second-order terms in the trailing wake boundary condition, 
and viscous effects modeled using integral boundary layer equations 14,15 solved 
simultaneously with the governing small perturbation potential flow equation with 
turbulent closure given by the dissipation integral relations of Drela. 16 The ASP theory 
was shown to be mathematically more appropriate and computationally more accurate 
than the classical TSP theories. Therefore, the ASP theory was used as the basis for a 
new computer program called ASP3D (Advanced Small Perturbation - 3D), which 
involves either AF1- or AF2-type approximate factorization algorithms and a FAS (full 
approximation scheme) multigrid procedure for the solution of the ASP potential flow 
equation and associated boundary conditions. The purpose of the program is for unsteady 
aerodynamic and aeroelastic analyses, especially in the nonlinear transonic flight regime. 
The ASP3D code can treat aircraft configurations involving multiple lifting surfaces with 
leading and trailing edge control surfaces and a fuselage. Steady-state aerodynamic 
results obtained using the ASP3D code were presented by Batina, thus demonstrating 
various options and the general utility of the new program. Comparisons were made with 
alternative calculations and experimental data to demonstrate the accuracy of the ASP 
approach for steady-state applications. The purpose of the present paper, though, is to 
present unsteady aerodynamic and aeroelastic applications obtained using the ASP3D 
code. Several cases were considered to assess the time-dependent capability and 
demonstrate additional features of ASP3D. Comparisons are made with alternative 
calculations and experimental data to evaluate the accuracy and utility of the code. 
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Overview of the ASP3D Program 
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The ASP3D code was developed over the last decade by the author during his 
employment as a Senior Research Scientist within the Configuration Aerodynamics 
Branch at NASA Langley Research Center. The development of the code resulted from a 
research study to determine the specific causes of the inaccuracies and inefficiencies 
inherent to CAP-TSD. That work revealed that the underlying small perturbation 
potential flow theory in CAP-TSD was deficient in several respects. An advanced small 
perturbation theory was developed to alleviate or eliminate the inaccuracies identified 
in CAP-TSD. The ASP theory was shown subsequently to be mathematically more 
appropriate and computationally more accurate than the classical TSP theories such as 
that embodied within CAP-TSD. Consequently, the ASP theory was used as the basis for 
the new ASP3D program, which in contrast with CAP-TSD involves modern 
computational fluid dynamics concepts such as a finite volume spatial discretization, 
sharper shock capture, dissipation integral boundary layer modeling, alternative AF2-type 
approximate factorization algorithm, and FAS (full approximation scheme) multigrid 
solution convergence acceleration. 

The ASP3D Computer Program 

The ASP3D computer program involves either AF1- or AF2-type approximate 
factorization algorithms and a FAS multigrid procedure for the solution of the advance 
small perturbation potential flow equation and associated boundary conditions. The 
ASP3D code can treat aircraft configurations involving multiple lifting surfaces with 
leading and trailing edge control surfaces, and a fuselage. The new code is described 
briefly here along with relevant equations. The interested reader is referred to References 
13 and 17 for further details of the program. 

Governing Equation 

The three-dimensional general small perturbation equation may be written in 
Cartesian coordinates as 

*L + *L + *2 + *l = 0 

dt dx dy dz 

where the Cartesian fluxes are defined as 

fo=~A(/> t -B(/) x 

/, =c+ D<p x + + F<p] + G(t; + 1 Ml 

fi = (i+«&+f <>;](>, +f$ 


3 



with the constants defined as 


A = Ml B = 2Ml C = 1 D = l-Ml 

E = -Uy+\)Ml F = -Uy+i)Ml G = Uy-3)Ml H = -{y-\)Ml 

2 6 2 

The ASP3D code though solves the governing equation written in computational 
coordinates as 

4 + i + 4 + |3 =0 

a d£ drj d£ 

where the computational fluxes are defined as 
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The various geometric quantities used in these equations are computed using the exact 
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formulas given by Batina. 

- Entropy and Vorticity Effects 

Shock-generated entropy effects are incorporated within small perturbation codes 
such as CAP-TSD by first using the Prandtl relation (shock jump condition) to 
determine the velocity downstream of the shock wave from the upstream and sonic 
velocities. The upstream and downstream velocities are then used in the Rankine- 
Hugoniot shock relation 18 to determine the change in entropy across the shock wave. 19 
The resulting change in entropy is subsequently used in a Clebsch formulation to 
determine the vorticity downstream of the shock wave. ' The vorticity modifies the 
calculation of the velocity field in the downstream region. This is a common procedure 
that has been used in various computer codes especially at the full-potential equation 
level. However, the approach as applied to the small perturbation equation, with any of 
the available definitions for the streamwise mass flux/i including that of Williams, 23 does 
not conserve mass. 
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In contrast, the approach developed for the ASP3D program conserves mass by 
using a ratio of the streamwise flux evaluated using the upstream and downstream 
velocities. Specifically, the downstream perturbation velocity is first computed using the 
Prandtl relation (shock jump relation) given by 


(&) 2 


0 + M,omA 

, + (^a)| 


where the subscripts 1 and 2 represent stations that are upstream and downstream of the 
shock wave, respectively. The change in entropy across the shock is then computed 
using 

As=(y- l){l - /, [(&), ]/ /, [(A ) 2 ]} 

For steady flows the entropy is held constant along gridlines downstream of shock waves. 
For unsteady flows the entropy is convected downstream using 

— As + — As = 0 
dt dx 

The fluxes downstream of the shock are subsequently modified according to 


1 / nonisentropic 
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which conserves mass across the shock wave by design. 

The ASP3D code uses a Clebsch formulation to compute the shock- generated 
vorticity. In brief, the streamwise velocities downstream of shocks are computed using 


Ox) 


rotational 


Ox) 


irrotational 


A 5 


Also, the second-order terms in the trailing wake boundary condition were found to be 

IV 

not insignificant, and thus, they were incorporated into the ASP3D solution procedure. 

- Viscous Effects 


The ASP3D viscous method involves solving the unsteady boundary layer equations 
simultaneously with the outer potential flow solution so that no interaction law coupling 
the inner and outer solutions is required. The combined solution procedure involves an 
implicit block tridiagonal inversion for all of the cells along the surface and trailing wake. 
Unlike the CAP-TSD viscous capability, 24 26 the ASP capability uses exact formulas for 
edge quantities and exact boundary conditions along surfaces and wakes. Smoothing of 
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edge quantities and limiters are not required for stability, and no arbitrary or free 
parameters are necessary to tune the procedure. 

Specifically in the ASP3D viscous capability, the exact edge quantities are used for 
the calculation of velocity u e , density p e , temperature T e , local Mach number M e , and the 
coefficient of viscosity p ( . defined by 


u e = 1 + fix 
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7 - 1 
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where S is Sutherland’s coefficient and T,„ is the freestream temperature. 

For the ASP3D capability, the time-dependent integral boundary layer (IBL) and lag 
equations may be written as a system of equations in the form 
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where the independent variables are the momentum thickness 6, the incompressible shape 
parameter Hk, and the square root of the shear stress coefficient c r . 

For solution, the integral boundary layer and lag equations are first premultiplied by 
[A]" 1 to diagonalize the time term as 
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The resulting IBL system is then linearized in a simple way and cast into the so-called 
delta-form for implicit solution as 
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where the viscous residual is defined as 
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The system of equations is directly coupled with the ASP equation for the outer 
potential flow for simultaneous implicit solution using a block tridiagonal matrix 
inversion procedure. This is done only for the cells adjacent to the upper and lower sides 
of the lifting surfaces and their wakes where the boundary layer is located. All other cells 
do not involve the boundary layer, and hence, only the ASP equation is solved in those 
cells. Dissipation integral relations are used for turbulent closure, similar to those 
reported by Drela. 16 

- Surface Boundary Condition 


In ASP3D the complete mass-consistent surface boundary condition including 
entropy, vorticity, and viscous effects is given by 


where 


and 
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The surface is represented by 

b(x,y,t)± &(x,t) — xa = 0 

where the ± is for the upper and lower surfaces, respectively, since the boundary layer 
displacement thickness 5 is defined as a positive quantity on both surfaces. 
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Note that the spanwise slopes b y are included in the complete surface boundary 
condition for three-dimensional applications, which may be important for swept wings 
near the leading edge and also in the wing tip region. Furthermore, the streamwise and 
spanwise slopes are computed internal to ASP3D from the ordinates of the wing 
geometry by simple second-order accurate central finite differences, effectively the same 
as is done in advanced codes with body fitted meshes. This eliminates the differentiation 
of the spline of the wing geometry normally performed as a preprocessing step to 
determine the airfoil slopes. It thus makes the input to the ASP3D code simpler than 
most other small perturbation codes such as CAP-TSD. 

- Wake Boundary Condition 

Viscous effects are included in the wake modeling in ASP3D through the boundary 
condition 

A <p : =±-(fS) + <S* 

8 v x 

across the wake where / and g were defined previously. The boundary condition is 
incorporated numerically using 


a) = AAAi±Q + l 
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Shock Capturing 

The AF1 and AF2 algorithms of the ASP3D code use one of three approaches to 
model supersonic regions of the flow including: (1) Godunov, (2) Engquist-Osher (E- 
O), and (3) Murman-Cole 29 (M-C) type-dependent switches. The first two switches 
satisfy an entropy inequality by design. This means that their use precludes the 
possibility of computing nonphysical expansion shocks that can occur with the Murman- 
Cole switch, especially in two-dimensional applications. With Godunov switching, 
shocks are captured sharply, usually with only one interior state. Therefore, the Godunov 
switch is preferable over the E-0 and M-C switches, and it is the default switch within 
the ASP3D code. 

Solution Algorithm 

The ASP3D code involves a modern cell-centered finite-volume discretization, 
which allows exact planform treatment including leading edges, wing tips, and control 
surface edges. The advanced small perturbation equation is solved numerically using 
either an AF1 approximate factorization algorithm or an AF2-type algorithm. The AF1 
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algorithm within the ASP3D code is the cell-centered finite-volume version of the AF1 

n 

algorithm developed for CAP-TSD, which is described by Batina. However, the AF1 

on 

algorithm is susceptible to the so-called moving shock instability if shock waves move 
more than one cell per time step. 

The AF2 algorithm within the ASP3D code is described in general by 


LgL/jL^At/) = -R(t/>) 


The general form of the AF2 algorithm appears to be similar to that of the AF1 algorithm 
but there are two fundamental differences. First, the operators are applied in the reverse 
order (here, vertical and spanwise sweeps are performed first, with the streamwise sweep 
performed last) and a streamwise implicit temporal damping term is included on the right 
hand side of the equation to eliminate the occurrence of the moving shock instability. 
The AF2 algorithm is defined specifically by 



S c g 3 S c 
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where the flux Jacobians are defined by Batina. Although not obvious, the temporal 
damping term is treated implicitly since it involves information from the second 
intermediate sweep but at the previous grid plane i - 1, and 


At 


a = — + 



2A t ) 


The parameter y/ is set equal to zero for steady applications and set equal to unity for 
unsteady applications. The parameter co is an over-relaxation factor. The AF2 algorithm 
involves both physical (At) and computational (A r) time steps. The computational step is 
selected for rapid convergence and it is usually assigned a large value. Hence the 
solution can be advanced using either a constant computational time step or a constant 
CFL number (variable time step). The physical time step is used only for unsteady 
applications, and it is selected based on the physics of the problem being considered. 

The AF2 approach is the preferred algorithm. The AF2 scheme is more robust than 
the AF1 scheme since it cannot fail due to the moving shock instability by design, 
because the temporal damping term eliminates the instability. Hence convergence to 
steady state is achieved with large time steps or large CFL numbers, independent of the 
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moving shock instability that plagues the AF1 algorithm. The AF2 scheme also may be 
used for unsteady calculations when local iteration is applied. Therefore, the physical 
time step may be selected based on the physics of the problem being solved, rather than 
on numerical stability considerations. In either steady or unsteady applications, the 
temporal damping term vanishes in the steady- state limit or within the convergence 
process of performing subiterations (for unsteady calculations). 

Multigrid Implementation 

The AF1 algorithm is not appropriate for multigrid implementation because the rapid 
rate of convergence of the multigrid procedure would result in the moving shock 
instability. The AF2 algorithm of the ASP3D code is amenable to multigrid 
implementation because the temporal damping term prevents the moving shock 
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instability from occurring. Therefore the multigrid procedure may be used for steady or 
iterative unsteady applications with ASP3D. The procedure is implemented in FAS (full 
approximation scheme) form, 32 which is applicable to nonlinear governing equations. As 
with any multigrid procedure, high-frequency errors are damped on the fine mesh and 
low-frequency errors are damped on the coarser meshes. Either V- or W-cycles may be 
selected and full multigrid is available. W-cycles have been found to be slightly more 
effective for solution convergence, although W-cycles are slightly more expensive than 
V-cycles. 33 

Specifically, the calculation on the finest mesh (N) smoothes high frequency errors 
using the AF2 algorithm written in general form 

Ln(0n) = 


Subsequent calculations on the next coarser mesh (AM) involve fine mesh residual 
injection 


^v-i(^v-i) - ^jv-i In (R n )+ L n _ { I n (0v) 


where the I functions are restriction operators for transferring the velocity potential and 
residual from the fine mesh to the coarse mesh. The procedure is generally repeated for 
additional coarse meshes (N-2, N- 3, . . .). 


The restriction operator for the velocity potential is a volume-weighted operator 
defined by 


0N-1 



ZV0N 

TV 


and the restriction operator for residuals is defined by 


I%~ l R N =TR 


N 


where the summations above are taken over all of the fine mesh cells that make up the 
coarse mesh cell. The interpolation operator for the velocity potential that is required to 


10 



transfer information back from coarse to fine mesh is a trilinear interpolation operator 
defined symbolically by 

0N 0N + IN-10N-1 


- Numerical Accuracy 

In the ASP3D code geometric quantities such as areas, volumes, and direction 
cosines are computed using exact formulas. The velocity components are computed 
using similar formulas as used in CAP-TSD, but in ASP3D, they are exact for a linear 
test function on a general grid. Generally speaking, the ASP3D algorithm is second- 
order- accurate within subsonic cells and first-order accurate within supersonic cells, 
although a second-order-accurate supersonic capability is available as a user-selected 
option. 

Also in CAP-TSD the (p xt term in the governing equation is computed using first- 
order- accurate backward differencing for stability reasons. However, for unsteady 
applications at higher reduced frequencies encountered with control surfaces or higher 
wave number phenomena such as Kutta waves, this differencing may be too dispersive. 
In contrast, first-, second-, and third-order-accurate backward differencing operators are 
available for the (j) x , term in ASP3D as user-selected options to improve the accuracy for 
such applications. 

- Structural Equations of Motion 

Similar to the CAP-TSD code, the flow equations of ASP3D are coupled with the 
structural equations of motion for simultaneous time- integration for aeroelastic analysis. 
Specifically, such analysis involves the coupling of the aerodynamics with the structural 
characteristics of the configuration under consideration. The resulting equations of 
motion for a time domain or time-marching aeroelastic solution are based on the aircraft 
natural vibration modes. These equations are integrated in time along with the finite 
volume solution of the flow field in ASP3D. Initial conditions for each structural mode 
are specified and free decay transients are computed. Initial velocities in one or more 
modes, rather than displacements, have been found to be superior in avoiding 
nonphysical, strictly numerical transients and their possibly associated instabilities. 
Aeroelastic stability is then deduced from the free decay records or time histories. This is 
a fairly standard procedure for aeroelastic analyses with the small perturbation computer 
codes, and thus, the interested reader is referred to Cunningham, et al. 2 or Bennett, et al. 3 
for further details including equations and representative calculations. 

Results and Discussion 

Results computed using the ASP3D computer program are presented for several 
unsteady aerodynamic and aeroelastic cases to assess the time-dependent capability and 
demonstrate various features of the code. The cases considered include: (1) the NACA 
0012 airfoil undergoing a forced harmonic pitching motion about the quarter chord at 
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transonic conditions, 34 (2) a thickening-thinning parabolic arc airfoil at a transonic Mach 
number, 35 (3) an investigation of wave propagation characteristics with emphasis on high 
wave number applications, 36 ' 37 (4) aeroelastic transient calculations for the NACA 0012 
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airfoil at values of the dynamic pressure below, near, and above the flutter value, ’ and 
(5) the F-5 fighter wing undergoing a rigid pitching motion about the wing root mid- 
chord axis. 40 All of the calculations are inviscid applications of the code, in part, because 
the viscous capability was not implemented at the time. 

Flow calculations for the two-dimensional cases were performed using a 193 x 129 
Cartesian finite volume mesh as shown in Figure 1. The total mesh is shown in Figure 
1(a) and a near field view is shown in Figure 1(b). The mesh had a twenty chord extent 
from the airfoil in both the x- and --directions with a spacing at the leading and trailing 
edges of Ax = 0.015c. There is only one cell in the spanwise direction for airfoil 
applications, which effectively produces a two-dimensional solution. Multigrid 
calculations were performed using this mesh (termed Mesh 5) and four subsets of the 
mesh (Meshes 1 through 4) determined by simply deleting every other grid line. The 
coarser meshes contain 97 x 65 (Mesh 4), 49 x 33 (Mesh 3), 25 x 17 (Mesh 2), and 13x9 
(Mesh 1) points in the x- and z-directions, respectively. This is done automatically 
within the ASP3D code and thus is not something that the user needs to do during mesh 
generation. Flow calculations for the F-5 fighter wing were performed using a 97 x 25 x 
65 finite volume mesh as discussed below. 

NACA 0012 Airfoil Harmonic Pitching Results 

The first case considered is that of the NACA 0012 airfoil at a mean angle of attack 
of oto = 0.016° and a freestream Mach number of M x = 0.755. 34 The airfoil was forced to 
pitch harmonically about the quarter chord with an amplitude of QC\ = 2.51° at a reduced 
frequency based on semichord of k = 0.0814. This is a well- studied nonlinear unsteady 
transonic case 34 that involves moving shock waves on both the upper and lower surfaces 
of the airfoil produced by the relatively large amplitude of motion. The shock waves 
move approximately twenty-five percent chord and periodically appear and disappear 
during a cycle of motion. The case was chosen to demonstrate the superior numerical 
stability of the iterative multigrid capability of the ASP3D code for unsteady 
applications, whereby the computational time step is selected to be very large for rapid 
convergence and the physical time step is selected based on the physics of the problem 
(temporal accuracy). In fact the physical step size can be selected to be very large as 
well, without creating numerical difficulties. 

To demonstrate the algorithm robustness, iterative multigrid calculations were 
performed for the NACA 0012 airfoil oscillating harmonically for four cycles of motion 
to ensure periodicity of the last cycle. The CFL number was selected to be 10 (which 
corresponds to a spatially varying computational step size) and the over-relaxation factor 
(co) was set equal to 1.6. The convergence criteria on a per time step basis was selected 
as two orders of magnitude, which produced an L 2 -norm of the residual of approximately 
10 at each time step. Calculations were performed first using only 32 steps per cycle of 
motion corresponding to a relatively large physical step size of At = 1.20608. This step 
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size corresponds to more than one chord of airfoil travel per time step. The forced 
pitching motion a(t) is shown in Figure 2(a) and the lift and moment coefficient 
responses, c\(t) and c m (t), respectively, are shown together in Figure 2(b). The 
coefficients indicate that the lift is a relatively linear response that slightly lags the 
pitching motion and the moment is a highly nonlinear response that lags the motion much 
more than the lift. Figures 2(c) and 2(d) show the iterative multigrid residual history and 
the number of supersonic points, respectively, during the last cycle of motion. The 
residual history essentially represents thirty-two successive multigrid calculations, one 
per physical time step, wherein the multigrid procedure is used to rapidly reduce the 
solution residual to the desired level, thus preserving temporal accuracy. By iterating on 
the solution at each time step, the linearization and factorization errors inherent in the 
AF2 algorithm are reduced to an acceptable level (~10 3 ) to produce a time accurate 
solution. The technique is analogous to performing subiterations per time step with CAP- 
TSD to ensure time accuracy, but the iterative multigrid procedure in ASP3D is cheaper 
computationally because the iterations are performed on coarser grids within the overall 
multigrid capability. The ASP3D calculation required about 180 iterations of the 
multigrid procedure, but since most of the work is done on the coarser Meshes 1 through 
4 using W-cycles, rather than on the finest Mesh 5, the total effort is less than iterating 
using only the finest mesh. The number of supersonic points during the last cycle of 
pitching motion (Figure 2(d)) is cyclic, as expected, with a larger value in the first half of 
the cycle than the value attained in the second half because of the small positive mean 
angle of attack (ocq = 0.016°). 

To further test and demonstrate the algorithm robustness of the ASP3D code, 
calculations for the pitching NACA 0012 airfoil were successively repeated using 
multiples of the original physical step size including At = 2.41216, 4.82432, and 9.64864, 
which correspond to 16, 8, and 4 steps per cycle (s/c) of motion, respectively. These 
values represent huge step sizes that are not usually numerically stable with other TSP 
computer codes. Admittedly, 16 and 8 steps per cycle barely describe the pitching 
motion in discrete form, and 4 steps per cycle of motion is the smallest value possible, 
which actually corresponds to a saw tooth pitching motion rather than a sinusoidal 
motion. Be that as it may, all of the calculations were numerically stable, and the results 
for 16, 8, and 4 s/c are presented in Figures 3, 4, and 5. The results are presented in the 
same format as the results of Figure 2 computed using 32 steps per cycle. Specifically, 
the pitching motion, the lift and moment coefficient time responses, the last cycle 
multigrid residual history, and the last cycle number of supersonic points are given in 
sub-figures (a), (b), (c), and (d), respectively. Comparisons of the results from the 16 s/c 
calculation (Figure 3) with the results obtained using 32 s/c (Figure 2) indicate that the 
lift response is very similar and the moment response is not as smooth, possibly because 
it is difficult to plot the nonlinear moment time history with only 16 straight line 
segments per cycle. The last cycle residual history shows very similar convergence to the 
32 s/c case but only 96 iterations of the multigrid procedure were required to obtain 
results for 16 s/c, which is approximately half the work required for the 32 s/c results. 

Comparisons of the results obtained using 8 s/c (Figures 4(a) - 4(d)) and 4 s/c 
(Figures 5(a) - 5(d)) with those using 32 s/c (Figures 2(a) - 2(d)) indicate a degradation 
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in the accuracy of the responses, especially the moment coefficient, but the calculations 
are numerically stable even though the corresponding step sizes are huge. The number of 
multigrid iterations required to obtain results with 8 s/c and 4 s/c are 48 and 24, 
respectively, as indicated in Figures 4(c) and 5(c). These values are exactly half and one- 
quarter the work required to obtain the results with 16 s/c (96 iterations), which indicates 
that there is a linear relationship between the number of steps per cycle and the number of 
multigrid iterations required to obtain results. For example, halving the number of steps 
per cycle in the calculation cuts the total work in half. 

Composite hysteresis plots of the responses in the last cycle of motion are plotted in 
Figure 6. Figure 6(a) shows the lift coefficient c\(t) versus the pitching motion aft) and 
Figure 6(b) shows the moment coefficient c m ( t) versus the pitching motion aft). The 
results from all four sets of calculations are presented corresponding to 32, 16, 8, and 4 
s/c of motion. The responses obtained using 32 and 16 s/c compare closely (the discrete 
points are connected with straight line segments) while the responses obtained using 8 
and 4 s/c are less accurate because the corresponding step sizes are too large. However, 
for easier cases involving smaller amplitudes of motion and less shock motion including 
flutter onset calculations (due to a small perturbation of the flow) the larger time steps 
may produce solutions of acceptable accuracy at a fraction of the cost. In any case, use 
of larger time steps does not create numerical stability problems, and therefore the analyst 
must make a prudent selection of step size balancing accuracy and cost of calculation. 

Thickening- Thinning Parabolic Arc Airfoil Results 

The second case considered is that of a parabolic arc airfoil with a time-varying 
midchord thickness at M„ <, = 0.85 and Oq = 0°. The temporal variation in thickness is 
defined mathematically and graphically in Figure 7. For example, beginning as a flat 
plate at t = 0, the thickness of the parabolic arc airfoil increases smoothly until it reaches 
maximum thickness of 10% (maximum thickness to chord ratio) at t = 15. After t = 15 
the airfoil thins until it reaches zero thickness at t = 30. This case is a simple two- 
dimensional helicopter rotor simulation that has been used by other researchers for code 
development purposes. 35 Similarly in the present study, the problem was solved using 
the ASP3D iterative multigrid procedure involving the five meshes discussed previously. 
The calculations were performed using a CFL number of 10 and a relatively large 
physical step size of At = 1.0 for a total of 32 steps. The step size of At = 1.0 corresponds 
to one chord of airfoil travel per time step. The convergence criteria on a per time step 
basis was selected as four orders of magnitude, which produced an L 2 -norm of the 
residual of approximately 10 5 at each time step. The resulting residual history and 
number of supersonic points are shown in Figures 8(a) and 8(b), respectively. The 
residual history (Figure 8(a)) indicates that 184 iterations were required to advance the 
solution 32 time steps and that each time step required about the same amount of work (5 
or 6 iterations) independent of whether the flow was subsonic or transonic. The number 
of supersonic points (Figure 8(b)) indicates that the solution goes from subsonic to 
transonic flow at time step 7, with a maximum of about 1250 supersonic points near time 
step 18 or 19, and back to subsonic flow at time step 27. 
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For comparison purposes the calculations were repeated with everything the same 
except for using only the finest Mesh 5. This calculation is referred to as an iterative 
unigrid calculation (no multigrid) to contrast it with the iterative multigrid calculation of 
Figure 8. The procedure then only iterates with the solution on the finest mesh, similar to 
that which is done in the CAP-TSD code when subiterations are applied on a per time 
step basis. The resulting residual history and number of supersonic points are presented 
in Figures 9(a) and 9(b), respectively. The residual history (Figure 9(a)) indicates that 
over 2400 unigrid iterations were required to achieve a similar convergence level as the 
iterative multigrid calculation (L 2 -norm of the residual of approximately 10 5 at each time 
step) and clearly more work (number of iterations) is required when the flow is transonic 
than when the flow is subsonic. Thus the iterative multigrid calculation (184 total 
iterations) is far more efficient than the iterative unigrid calculation (over 2400 total 
iterations), because most of the iterative multigrid work is done on coarser meshes 
(Meshes 1 through 4). The number of supersonic points (Figure 9(b)) is similar to that 
shown from the multigrid calculation (Figure 8(b)), such as attaining a maximum of 
about 1250 supersonic points, since the two flow field solutions should be almost 
identical. 

Instantaneous solutions from the iterative multigrid calculations at t = 8.5, 11.5, 
18.25, 26.875, 29.0, and 32.0 are shown in Figures 10, 11, 12, 13, 14, and 15, 
respectively. These six points in time were selected for later comparison with the 
published results of Goorjian. 35 The top pail (a) of each figure shows the surface 
pressure coefficient distribution and the bottom part (b) of each figure shows the near 
field pressure coefficient contour lines. For t = 29.0 and 32.0, the results were taken 
directly from the multigrid calculations discussed above since the step size was A t = 1.0 
(exactly 29 and 32 total time steps, respectively). To avoid solution interpolation to 
obtain the results for t = 8.5, 11.5, 18.25, and 26.875, the calculations were repeated four 
times with step sizes equal to At = 8.5/8, 11.5/11, 18.25/18, and 26.875/26, respectively, 
using 8, 11, 18, and 26 total time steps. 

At t = 0 the airfoil has zero thickness (flat plate) so the undisturbed flow is the 
correct starting solution for the thickening-thinning airfoil problem. As the thickness of 
the parabolic arc airfoil increases in time the solution develops such that, for example at t 
= 8.5 (Figures 10(a) and 10(b)), the flow is relatively smooth with very small embedded 
supersonic regions forming on each side of the airfoil with weak shock waves beginning 
to develop near 76% chord. As the airfoil further thickens to t = 11.5 (Figures 11(a) and 
11(b)) the shock waves grow to moderate strength near 83% chord. The airfoil reaches 
its maximum thickness of 10% at t = 15, so by t = 18.25 (Figures 12(a) and 12(b)), the 
embedded supersonic regions are approximately at their largest (1250 supersonic points) 
and the relatively strong shock waves are located near 91% chord. As the airfoil becomes 
thinner in time the flow becomes subsonic again and the upper and lower surface shock 
waves become upstream propagating pressure waves near midchord as indicated at t = 
26.875 (Figures 13(a) and 13(b)). At t = 29.0 (Figures 14(a) and 14(b)) the pressure 
waves have traveled upstream to roughly the airfoil quarter-chord location and by t = 
32.0 (Figures 15(a) and 15(b)), when the airfoil has zero thickness again, the pressure 
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waves on the upper and lower surfaces have moved forward of the airfoil leading edge 
thus forming a single wave. 

Comparisons of surface pressure coefficients for the thickening-thinning parabolic 
arc airfoil case are presented in Figure 16. The ASP3D pressure distributions at t = 8.5, 
11.5, and 18.25 are shown in Figure 16(a) with corresponding full-potential (FP) and 
TSD pressure distributions from Goorjian shown in Figure 16(b). The FP and TSD 
results of Reference 35 were obtained using step sizes of 0.0625 and 0.125, respectively, 
whereas the ASP3D results were obtained using step sizes equal to or greater than 1.0. 
The ASP3D pressures (Figure 16(a)) are in good agreement with the pressures reported 

ic 

by Goorjjian, (Figure 16(b)) especially in comparison with the full-potential pressure 
results. Specifically at t = 18.25 where the shock waves are strongest, the ASP3D shock 
location compares better with the FP result than the TSD result. Further pressure 
comparisons for t = 26.875, 29.0, and 32.0 are presented in Figure 16(c) from ASP3D 
and in Figure 16(d) from Reference 35. The ASP3D pressures again agree best with the 
full-potential pressure distributions, although all of the results indicate that the pressure 
waves are smeared as they propagate upstream off the airfoil leading edge. 

- Linear Wave Propagation Results 

To understand the importance of accurate wave propagation for unsteady 
aerodynamic flows, it is instructive to relate some of the work reported over thirty years 
ago by Tijdeman and Zwaan 36 at the NLR in Amsterdam. Reference 36 contains a 
discussion of the propagation of pressure waves in high subsonic and transonic flows. 
The authors looked at disturbance wave fronts generated by a harmonically pulsating 
source at the three-quarter-chord location (control surface hinge line) along an airfoil, 
presented here in Figure 17 taken from their report. Figure 17(a) shows upstream 
propagating wave fronts for a subsonic freestream Mach number of A7,„ = 0.8 and Figure 
17(b) similarly shows wave fronts for a transonic Mach number of M „ „ = 0.875. The 
latter flow involves an embedded supersonic region on the airfoil terminated by a shock 
wave near midchord. In the upper parts of the two figures the positions of the wave 
fronts are drawn at equal time intervals for the slowly moving receding waves. The 
lower parts of the figures show the same wave fronts, but for a uniform flow (i.e., without 
the airfoil) as a reference. At the subsonic freestream Mach number of = 0.8 (Figure 
17(a)), the actual (upper) wave fronts move a little slower than the uniform flow (lower) 
as suggested by the closer spacing of the curves representing the fronts. The spacing of 
the wave fronts is a measure of the intensity of the local pressure gradients. Also the 
wave fronts incline forward slightly due to the vertical velocity gradient. At the transonic 
Mach number of M „ = 0.875 (Figure 17(b)), the lower part of the actual wave fronts 
coalesce into the shock wave. The upper parts of the wave fronts rotate about the top of 
the shock and enter the supersonic region at an oblique angle, inclined in comparison to 
the uniform flow. The waves are slowed significantly near the shock wave thus 
producing a strong pressure gradient there, as expected. Also, the longer time required 
for the wave fronts to travel around the shock wave leads to a significant time lag in the 
airfoil aerodynamic response. Although not demonstrated here, this time lag is generally 
regarded as the physical mechanism responsible for the so-called transonic flutter dip. 
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Tijdeman and Zwaan 36 also presented wind tunnel results in the form of 
shadowgraph photographs of high-speed flows about an airfoil reproduced here in Figure 
18. The upstream propagating wave fronts are generated by vortex streets emanating 
from the aft end of the airfoil with its trailing edge control surface removed. Figure 18(a) 
shows the shadowgraph for the subsonic Mach number of M <*, = 0.8 wherein the receding 
wave fronts have relatively uniform spacing and they are slightly inclined as in the 
diagram of Figure 17(a). Figure 18(b) shows the shadowgraph for the transonic Mach 
number of M„ „ = 0.875 wherein the lower parts of the receding wave fronts coalesce near 
the foot of the shock wave and the upper parts of the fronts rotate around the top of the 
shock to enter the supersonic region, similar to that diagrammed in Figure 17(b). These 
figures illustrate some of the flow field complexities for unsteady transonic aerodynamics 
that computational methods must be able to model accurately. 

It also is instructive to examine the simpler linear theory wave propagation 
characteristics of the upstream and downstream propagating wave fronts in a uniform 
flow and the implications for resolving such waves using discrete computational 
methods. For example, Figures 19(a) and 19(b) show wave numbers for upstream and 
downstream propagating wave fronts, respectively, as a function of Mach number for 
various values of reduced frequency k. The wave number associated with upstream 
propagating wave fronts is defined by 


wave number = 


Af 


1 ~M„ 


k 


and the wave number associated with downstream propagating wave fronts is defined by 

wave number = ^°° — k 
l+M„ 

According to Rowe and Cunningham, “The wave number is a useful parameter for 
assessing the waviness of pressure distributions and estimating relative difficulties in 
predicting converged values of generalized forces.” Computational fluid dynamics 
(CFD) methods such as CAP-TSD or ASP3D are typically capable of computing 
propagating wave fronts for wave numbers less than about 5.0 depending upon the 
density of the mesh and the size of the physical time step. Therefore the computation of 
most downstream propagating waves (shown in Figure 19(b)) with associated wave 
numbers less than about 2.0 is not difficult because the waves are moving very fast (with 
the freestream), although relatively small time steps may be required to accurately 
describe downstream running waves at higher reduced frequency. In contrast, the 
computation of upstream propagating wave fronts is difficult at higher reduced 
frequencies because the associated wave number is much higher, as shown in Figure 
19(a). Pressure distributions with wave numbers higher than about 5.0 have significant 
waviness that can be difficult to capture with CFD methods as discussed below. As the 
Mach number increases toward unity from the subsonic side the problem is exacerbated, 
because the waves are moving increasingly slowly (against the freestream). The wave 
number curves in Figure 19(a) asymptote to infinity at Mach one because the wave fronts 
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there are theoretically standing waves. And for supersonic freestream there can be no 
upstream propagation. 

Another way of looking at the linear theory wave propagation characteristics is by 
plotting the reduced frequency of the upstream and downstream propagating wave fronts 
as functions of Mach number for various values of the wavelength A as presented in 
Figures 20(a) and 20(b), respectively. The upstream values are computed by 

k = l-M„ n 

X 


and the downstream values are computed by 

1 + k 

A = v 

X 

For the upstream propagating wave fronts (Figure 20(a)) the reduced frequency and 
wavelength are inversely proportional. As the Mach number approaches unity from the 
subsonic side the reduced frequency is reduced for a given wavelength until the singular 
limit at Mach one is reached where a standing wave can occur. The small wavelength 
fronts that may be regarded as being outside the range of practical reduced frequencies 
for Mach numbers of, say, = 0.8, are in fact within the range encountered in 
aeroelasticity (approximately 0 < k < 2) at higher Mach numbers of 0.9 or 0.95. 
Therefore the accurate resolution of these small wavelength wave fronts should be of 
concern for CFD development. Of course for supersonic freestream there is no upstream 
wave propagation. For the downstream propagating wave fronts (Figure 20(b)) the 
reduced frequency is again inversely proportional to the wavelength. But here, 
wavelengths less than about 4.0 result in reduced frequencies that are outside the typical 
range encountered in aeroelasticity. Conversely, wavelengths greater than 4.0 correspond 
to reduced frequencies less than 2.0 for any Mach number, but since these wavelengths 
are so large, they can be represented easily with most CFD codes given a sufficient grid 
and appropriate step size. Therefore it is the upstream propagating waves where the 
performance of these codes needs to be tested. 

Rowe and Cunningham presented lifting pressure coefficient distributions for 
unsteady cases involving high wave numbers. Similarly, Figure 21 shows the effects of 
Mach number on the lifting pressure coefficient distributions for a flat plate airfoil 
oscillating in pitch about the leading edge at k = 0.9. Figure 21(a) shows the real or in- 
phase component of the pressure and Figure 21(b) shows the imaginary or out-of-phase 
component of the pressure. The results were computed using the two-dimensional 
assumed-pressure-mode kernel-function linear-aerodynamics program of Bland. 41 The 
figures clearly show that as the freestream Mach number increases from 0.75 to 0.95 
more waves occur in the pressure distribution. For the Mach numbers of 0.75, 0.80, 0.85, 
0.90, and 0.95, the corresponding upstream propagating wave numbers are 2.7, 3.6, 5.1, 
8.1, and 19.0, respectively. 
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Figure 22(a) shows the lifting pressure coefficient distribution computed using the 
ASP3D code for the flat plate airfoil oscillating in pitch about the leading edge at k = 0.9. 
The freestream Mach number was selected as = 0.9 and the linear theory results from 
the previous Figure 21 that had a relatively high wave number of 8.1 are repeated here for 
comparison purposes. It also is a case studied in Reference 37. Figure 22(a) shows the 
real and imaginary ASP3D pressures, which agree with the linear theory pressures near 
the trailing edge, but do not agree with linear theory as the waves propagate upstream 
toward the leading edge. The ASP3D results were obtained using a first-order-accurate 
backward difference formula to compute the <p xt term. A backward difference formula is 
typically used to compute the (p xt term in TSP codes for stability reasons, and similar 
erroneous results were obtained using CAP-TSD. However in ASP3D, second- and 
third-order-accurate formulas have been successfully implemented to improve the 
temporal accuracy of the results for higher wave number cases. The pressure distribution 
for this case computed using the second-order-accurate formula is presented in Figure 
22(b). These results show a marked improvement in accuracy in comparison with the 
linear theory result. Results using the third-order-accurate formula are identical to 
plotting accuracy with the second-order-accurate results, and thus, they are not presented 
here. 


NACA 0012 Airfoil Aeroelastic Transient Results 

To assess and demonstrate the aeroelasticity capability in ASP3D, calculations were 
performed for the NACA 0012 airfoil at = 0.8 and Ofo = 0°, a case considered by 
Robinson, et al. using the CFL3D code run in Euler mode. Steady-state calculations 
were performed first, using ASP3D including entropy and vorticity effects, to provide the 
initial conditions for aeroelastic computation upon restart. The resulting steady pressure 
coefficient distribution is presented in Figure 23(a) with the corresponding CFL3D 
(Euler) pressure distribution from Ref. 38 presented in Figure 23(b). The ASP3D 
pressure distribution is in very good agreement with the CFL3D pressure distribution 
whereby the shock waves on the upper and lower surfaces of the airfoil are predicted to 
be near midchord. In the ASP3D calculation the shock waves are captured sharply 
because of the use of the Godunov type-dependent switch and the pressure levels are well 
predicted including those near the airfoil leading edge (stagnation). 

For aeroelastic analysis, the NACA 0012 airfoil was modeled structurally with two 
degrees of freedom, bending and torsion. The case considered is the much- studied Case 
“A” of Isogai which has normal modes similar to those of a streamwise section of a 
sweptback wing. The wind-off bending and torsion natural frequencies are 71.33 and 
535.65 rad/s, respectively. The pivot point for the bending mode is located 1.44 
chordlengths upstream of the leading edge of the airfoil. The pivot point for the torsion 
mode is 0.068 chordlengths forward of midchord. These mode shapes and natural 
frequencies were determined by performing a free vibration analysis with the aeroelastic 
equations written in the traditional form of plunge and pitch degrees of freedom. In this 
analysis, the following structural parameter values were used: a = -2.0, x a = 1.8, r a = 
1.865, = 100 rad/s, and = 100 rad/s. Also, the airfoil mass ratio was fl = 60. 

Generalized displacements corresponding to the bending and torsion modes are defined 
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as q\(t) and cj2(t), respectively. Initial conditions for the time-marching aeroelastic 
calculations were dq\(0)/dt = 2.0 and dcj2(0)/dt = 1.0. 

Aeroelastic calculations were performed for the NACA 0012 airfoil at M„ „ = 0.8 and 
O(o = 0° for several values of the nondimensional dynamic pressure including Q = 0.2, 
0.5, and 0.8, where Q = (U. <Jb cotdl jd ) 2 . In these calculations a smaller value of the step 
size was selected than in previous unsteady calculations, At = 0.1, to ensure that the 
higher frequency torsion mode would be resolved in the aeroelastic integration. Also the 
multigrid procedure was not exercised because the computational benefit would not be as 
large for such a step size. The generalized displacements for Q = 0.2, 0.5, and 0.8 are 
presented in Figures 24, 25, and 26, respectively. For Q = 0.2 (Figure 24), the 
generalized displacements indicate that the aeroelastic system is stable because the 
transients are converging in time. The qi(t) generalized displacement is dominated by the 
lower frequency bending mode and the cj2(t) generalized displacement contains both 
modes including the higher frequency torsion mode. Use of the step size of At = 0.1 
resulted in approximately 133 s/c for the bending mode and only 20 s/c for the higher 
frequency torsion mode. Use of much larger time steps, although numerically stable, 
produces results that inaccurately describe the torsion mode. 

For Q = 0.5 (Figure 25), the generalized displacements indicate that the aeroelastic 
system is near-neutrally stable (flutter) since the dominant bending mode in both q\(l) 
and q2(t) is not damped. This correlates well with the reported results from Ref. 38 
wherein the transients computed using CFL3D (Euler) at Q = 0.5 were described as near- 
neutrally stable corresponding to flutter. Also the q2(t) generalized displacement 
indicates that the higher frequency torsion mode is more damped than at Q = 0.2. Finally 
for Q = 0.8 (Figure 26), the generalized displacements q\(t) and q 2 (t) are both diverging 
corresponding to an unstable aeroelastic system and the higher frequency torsion mode 
contained within q2(t) has even more damping than the Q = 0.2 and 0.5 transients. 

F-5 Fighter Wing Harmonic Pitching Results 

The F-5 fighter wing 40 has a leading edge sweep of 31.92 degrees, an aspect ratio of 
1.578, and a taper ratio of 0.283. The wing has a NACA 65A004.8 airfoil section that 
has a drooped leading edge. The wing was modeled for ASP3D applications using a 97 x 
25 x 65 finite volume mesh as shown in Figure 27. A near field view of the 97 x 25 
planform mesh is shown in Figure 27(a) and a near field view of the 97 x 65 root 
sectional mesh is shown in Figure 27(b). The meshes were generated using a program 43 
that was created to construct meshes specifically for use by the ASP3D code, and as such, 
it can generate meshes for aircraft configurations with multiple lifting surfaces, with 
multiple leading and trailing edge controls, and a fuselage. The mesh generator 
automatically creates the input file for ASP3D by using an extension of the polynomial 
mesh generation procedure developed by Bland. 44 Here though, the equations were 
generalized to the physical domain rather than the computational domain, which allows 
the creation of three-dimensional meshes about swept tapered planforms with the ability 
to control the spacing along the leading and trailing edges. This is especially an issue in 
the outboard region of the wing near the tip. For example, Figures 28(a) and 28(b) show 
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near field views of the planform mesh near the wing tip leading and trailing edges, 
respectively, for the F-5 wing. The figures show that the mesh was generated with a 
uniform spacing about the leading and tailing edges in the streamwise direction, and in 
fact, the spacing is also uniform across the tip in the spanwise direction. The planform 
mesh also emphasizes that all edges of the lifting surface are modeled with grid lines for 
ASP3D application. This is in contrast with meshes generated for CAP-TSD application, 
which cannot have grid lines on leading edges or the edge of the wake (wing tip) because 
of singularities that occur along these lines. Details of the mesh generation procedures 
developed for the ASP3D computer program are presented by Batina. 43 

To demonstrate the multigrid procedures for a three dimensional case, calculations 
were performed using the 97 x 25 planform mesh (termed Mesh 4) and three subsets of 
the mesh (Meshes 1 through 3) determined by simply deleting every other grid line. This 
is done automatically within the ASP3D code and thus is not something that the user 
needs to do during mesh generation. The four planform meshes are shown in Figures 
29(a) to 29(d) including the total number of points in each direction for each mesh. The 
smallest mesh, for example, is a very coarse 13 x 4 mesh, which has only a few cells on 
the surface of the wing. Similarly, near field views of the corresponding four sectional 
meshes at the root of the F-5 wing are shown in Figures 30(a) through 30(d) with the 
airfoil slit in the center of each plot. The individual figures illustrate the relative mesh 
density near the airfoil surface along the wing root. 

To demonstrate the iterative multigrid capability for time-dependent applications in 
three dimensions, unsteady calculations were performed for the F-5 fighter wing pitching 
about the root midchord with an amplitude of a\ = 0.109° and a reduced frequency based 
on root semichord of k = 0.137. In these calculations the freestream Mach number was 
Moo = 0.899 and the mean angle of attack was (Xq = -0.005°. The wing was forced to 
oscillate harmonically in pitch for four cycles of motion as shown in Figure 31(a), to 
ensure periodicity of the flow field, with the resulting lift and moment coefficients shown 
in Figure 31(b). The calculations were performed using only sixteen steps per cycle of 
motion, corresponding to a very large time step of At = 1.433, wherein the multigrid 
procedure was used at each time step to reduce the solution residual by two orders of 
magnitude as shown in Figure 31(c). By iterating on the solution at each time step, the 
linearization and factorization errors inherent in the AF2 algorithm are reduced to an 
acceptable level (~10 4 ) to produce a time accurate solution. The technique is analogous 
to performing subiterations per time step with CAP-TSD to ensure time accuracy, but the 
iterative multigrid procedure in ASP3D is cheaper computationally because the iteration 
is performed on coarser grids within the overall multigrid capability. The corresponding 
number of supersonic points is presented in Figure 31(d). 

The resulting unsteady pressure coefficient distributions, normalized by the 
amplitude of motion and interpolated to the span stations where the NLR experimental 
pressure data 40 were measured, are presented in Figure 32. The unsteady pressure 
distributions on the upper surface of the wing are shown in Figure 32(a) and the unsteady 
pressure distributions on the lower surface of the wing are shown in Figure 32(b). The 
ASP3D pressure distributions were computed from the fourth cycle of motion and the 
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real and imaginary parts of the pressures are the in-phase and out-of-phase components of 
the pressure time histories, respectively. The unsteady pressures on the upper surface of 
the wing agree reasonably well with the experimental data 40 except near the shock wave. 
Here the magnitude of the pressure pulse due to the movement of the shock wave is 
slightly over predicted in the ASP3D calculation because of the neglect of viscous 
effects. On the lower surface of the wing (Figure 32(b)) the unsteady pressures are in 
very good agreement with the experimental data, with relatively uniform agreement with 
the data from the wing root to the wing tip. 

Conclusions 

A new computer program has been developed called ASP3D (Advanced Small 
Perturbation - 3D), which solves the small perturbation potential flow equation in an 
advanced form including mass-consistent surface and trailing wake boundary conditions, 
and entropy, vorticity, and viscous effects. The purpose of the program is for unsteady 
aerodynamic and aeroelastic analyses, especially in the nonlinear transonic flight regime. 
The program exploits the simplicity of stationary Cartesian meshes with the movement or 
deformation of the configuration under consideration incorporated into the solution 
algorithm through a planar surface boundary condition. The ASP3D code is the result of 
a decade of developmental work on improvements to the small perturbation formulation, 
performed while the author was employed as a Senior Research Scientist in the 
Configuration Aerodynamics Branch at the NASA Langley Research Center. The 
ASP3D code is a significant improvement to the state-of-the-art for transonic aeroelastic 
analyses over the CAP-TSD code (Computational Aeroelasticity Program - Transonic 
Small Disturbance), which was developed principally by the author in the mid-1980s. 

The ASP3D code involves a modern cell-centered finite-volume discretization, 
which allows exact planform treatment including leading edges, wing tips, and control 
surface edges. It can treat aircraft configurations involving multiple lifting surfaces with 
leading and trailing edge control surfaces, and a fuselage. The advanced-small- 
perturbation potential flow equation is solved numerically using either an AF1- or an 
AF2-type approximate factorization algorithm. The AF1 algorithm within the ASP3D 
code is the cell-centered finite-volume version of the AF1 finite-difference algorithm 
developed for CAP-TSD. However, the AF1 algorithm is susceptible to the so-called 
moving shock instability if shock waves move more than one grid point or cell per time 
step. In contrast, the AF2 scheme has a temporal damping term that eliminates the 
moving shock instability by design, and thus, the AF2 algorithm is more robust than the 
AF1 scheme. Thus, with the AF2 approach, convergence to steady state is achieved with 
large time steps or large CFL numbers. And because the moving shock instability does 
not occur, the AF2 algorithm of the ASP3D code is amenable to multigrid 
implementation. Therefore, a multigrid procedure was developed within ASP3D that 
may be used for steady or iterative unsteady applications. The procedure was 
implemented in FAS (full approximation scheme) form, which is applicable to nonlinear 
governing equations. The multigrid procedure may be used for unsteady calculations 
when local iteration is applied. Therefore, the physical time step may be selected for the 
accuracy of the problem under consideration independent of numerical stability issues, 
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and a large computational step size or CFL number may be selected for rapid 
convergence of the iterative procedure. 

Unsteady aerodynamic and aeroelastic applications of ASP3D were presented to 
assess the time dependent capability and demonstrate various features of the code. The 
cases considered included: (1) the NACA 0012 airfoil undergoing a forced harmonic 
pitching motion about the quarter chord at transonic conditions, (2) a thickening-thinning 
parabolic arc airfoil at a transonic Mach number, (3) an investigation of wave 
propagation characteristics with emphasis on high wave number applications, (4) 
aeroelastic transient calculations for the NACA 0012 airfoil at values of the dynamic 
pressure below, near, and above the flutter value, and (5) the F-5 fighter wing undergoing 
a rigid pitching motion about the wing root midchord axis. The results compared well 
with alternative methods and experimental data, and thus demonstrated the efficiency and 
utility of the ASP3D code for various unsteady aerodynamic and aeroelastic applications. 
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(a) Total mesh. 



x/c 


(b) Near field view. 

Figure 1 - Cartesian finite- volume mesh (193 x 129) with twenty-chord extent for 
airfoil unsteady aerodynamic and aeroelastic applications with the ASP3D 

computer program. 
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Figure 2 - Iterative multigrid calculation using 32 steps per cycle of motion for the 
NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 
Moo = 0.755, oq = 0.016°, a x = 2.51°, and k = 0.0814. 
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(d) Last cycle of the number of supersonic points. 


Figure 2 - Concluded. 
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(a) Forced harmonic pitching for four cycles of motion. 




Figure 3 - Iterative multigrid calculation using 16 steps per cycle of motion for the 
NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 
Moo = 0.755, oq = 0.016°, a x = 2.51°, and k = 0.0814. 


30 







16 steps/cycle 



Figure 3 - Concluded. 
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(a) Forced harmonic pitching for four cycles of motion. 




Figure 4 - Iterative multigrid calculation using 8 steps per cycle of motion for the 
NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 
Moo = 0.755, ct> = 0.016°, a, = 2.51°, and k = 0.0814. 
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Figure 4 - Concluded. 
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(a) Forced harmonic pitching for four cycles of motion. 




Figure 5 - Iterative multigrid calculation using 4 steps per cycle of motion for the 
NACA 0012 airfoil oscillating harmonically in pitch about the quarter chord at 
Moo = 0.755, ct> = 0.016°, a, = 2.51°, and k = 0.0814. 
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4 steps/cycle 



(d) Last cycle of the number of supersonic points. 


Figure 5 - Concluded. 
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(a) Lift coefficient. 



«(t) 


(b) Moment coefficient. 


Figure 6 - Effects of step size on the lift and moment coefficient responses in the last 
(fourth) cycle of motion for the NACA 0012 airfoil oscillating harmonically in pitch 
about the quarter chord at = 0.755, = 0.016°, a\ = 2.51°, and k = 0.0814. 
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Figure 7 - Definition of midchord thickness as a function of time for the thickening- 
thinning parabolic arc airfoil. 
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supersonic 



Figure 8 - Iterative multigrid calculation for the thickening-thinning parabolic arc 
airfoil problem at = 0.85 (32 total time steps). 
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(a) Residual history. 



iterations 

(b) Number of supersonic points. 

Figure 9 - Iterative unigrid calculation for the thickening-thinning parabolic arc 
airfoil problem at = 0.85 (32 total time steps). 
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(a) Pressure coefficient distribution. 



Figure 10 - Instantaneous results at t = 8.5 from the iterative multigrid 
calculation of the thickening-thinning parabolic arc airfoil at = 0.85. 
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(b) Near field pressure coefficient contour lines. 


Figure 11 - Instantaneous results at t = 11.5 from the iterative multigrid 
calculation of the thickening-thinning parabolic arc airfoil at M„ „ = 0.85. 
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(a) Pressure coefficient distribution. 
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(b) Near field pressure coefficient contour lines. 


Figure 12 - Instantaneous results at t = 18.25 from the iterative multigrid 
calculation of the thickening-thinning parabolic arc airfoil at = 0.85. 
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(a) Pressure coefficient distribution. 



x/c 

(b) Near field pressure coefficient contour lines. 


Figure 13 - Instantaneous results at t = 26.875 from the iterative multigrid 
calculation of the thickening-thinning parabolic arc airfoil at M„ „ = 0.85. 
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(a) Pressure coefficient distribution. 
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(b) Near field pressure coefficient contour lines. 

Figure 14 - Instantaneous results at t = 29.0 from the iterative multigrid 
calculation of the thickening-thinning parabolic arc airfoil at M„ „ = 0.85. 


44 





- 0.8 
- 0.6 
-0.4 
- 0.2 
0.0 
0.2 
0.4 
0.6 

-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 

x/c 

(a) Pressure coefficient distribution. 
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(b) Near field pressure coefficient contour lines. 


Figure 15 - Instantaneous results at t = 32.0 from the iterative multigrid 
calculation of the thickening-thinning parabolic arc airfoil at M„ „ = 0.85. 
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(a) ASP3D pressure coefficient distributions at t = 8.5, 11.5, and 18.25. 

FULL-POTENTIAL EQUATION 

Lt = 0.0625 

LOW-FREQUENCY, SMALL- 

DISTURBANCE, TRANSONIC 



(b) Corresponding FP and TSD pressure coefficients reported by Goorjian. 35 

Figure 16 - Pressure coefficient comparisons for the thickening-thinning parabolic 

arc airfoil at = 0.85. 
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(c) ASP3D pressure coefficient distributions at t = 26.875, 29.0, and 32.0 



X .5 


(d) Corresponding FP and TSD pressure coefficients reported by Goorjian. 3 


Figure 16 - Concluded. 




(a) Subsonic freestream Mach number of M «, = 0.8. 



(b) Transonic freestream Mach number of M„ 0 = 0.875. 

Figure 17 - NLR diagrams 36 of upstream propagation of disturbance wave fronts 
generated by a source at the control surface hinge line. 
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(a) Upstream propagation. 
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(b) Downstream propagation. 
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Figure 19 - Linear theory wave propagation characteristics as a function of Mach 
number for a range of values of reduced frequency k. 
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(a) Upstream propagation. 
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X = 2.0 
X = 1.0 


(b) Downstream propagation. 


Figure 20 - Linear theory wave propagation characteristics as a function of Mach 
number for various values of wavelength X. 
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(a) Real part. 
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(b) Imaginary part. 

Figure 21 - Effects of Mach number on linear theory lifting pressure coefficient 
distributions for a flat plate airfoil oscillating about the leading edge at k = 0.9. 
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(b) Second-order-accurate (j) xt term. 

Figure 22 - Effects of order of accuracy on ASP3D lifting pressure coefficient 
distributions for a flat plate airfoil oscillating about the 1. e. at k = 0.9 and = 0.9. 


53 






x/c 

(a) ASP3D pressure distribution including entropy and vorticity effects. 



(b) Euler (CFL3D) pressure distribution computed by Robinson, et al. 38 

Figure 23 - Steady pressure coefficient distributions for the NACA 0012 airfoil at 

Moo = 0.8 and a« = 0°. 
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Figure 24 - Generalized displacements corresponding to stable (converging) 
aeroelastic motion for the NACA 0012 airfoil at = 0.8, ceb = 0°, and Q = 0.2 for 
the structural parameter values of Isogai Case “A”. 
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Figure 25 - Generalized displacements corresponding to near-neutrally stable 
(flutter) aeroelastic motion for the NACA 0012 airfoil at = 0.8, Ofa = 0°, and 
Q = 0.5 for the structural parameter values of Isogai Case “A”. 
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Figure 26 - Generalized displacements corresponding to unstable (diverging) 
aeroelastic motion for the NACA 0012 airfoil at = 0.8, ab = 0°, and Q = 0.8 for 
the structural parameter values of Isogai Case “A”. 
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(b) Near field view of 97 x 65 root sectional mesh. 
Figure 27 - Finite volume meshes for the F-5 fighter wing. 
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(a) Mesh 4 (97x25). 


(b) Mesh 3 (49x13). 



(c) Mesh 2 (25x7). 


(d) Mesh 1 (13x4). 


Figure 29 - Near field view of finite volume planform meshes for the F-5 fighter 
wing for ASP3D multigrid calculations. 
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(a) Mesh 4 (97x65). 
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(b) Mesh 3 (49x33). 



(c) Mesh 2 (25 x 17). (d) Mesh 1 (13 x 9). 


Figure 30 - Near field view of finite volume root sectional meshes for the F-5 fighter 

wing for ASP3D multigrid calculations. 


61 



0.14 
0.07 
a(f) 0.00 
-0.07 
-0.14 

0 20 40 60 80 100 

time 

(a) Forced harmonic pitching for four cycles of motion. 
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(b) Lift and moment coefficient responses. 

Figure 31 - Iterative multigrid calculation using 16 steps per cycle of motion for the 
F-5 fighter wing pitching about the root midchord at = 0.899, 

(%) = -0.005°, a\ - 0.109°, and k = 0.137. 
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(c) Iterative multigrid residual history. 
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ASP3D inviscid solution 
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(a) Upper surface. 







(b) Lower surface. 


Figure 32 - ASP3D inviscid unsteady pressure coefficient comparisons with NLR 
experimental data 40 for the F-5 fighter wing pitching at M„ „ = 0.899, 

= -0.005°, a\ = 0.109°, and k = 0.137. 
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